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We study statistical mechanics of the self-gravitating system applying the cluster expansion method 
developed in solid state physics. By summing infinite series of diagrams, we derive a complex free 
energy whose imaginary part is related to the relaxation time of the system. Summation of another 
series yields two-point correlation function whose correlation length is essentially given by the Jeans 
wavelength of the system. 
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a. Introduction Variety of structures found in the 
present Universe are thought to be seeded by quan¬ 
tum fluctuations in the very early Universe and evolved 
through the gravitational interaction in the expanding 
Universe aided by yet unknown Dark Matter. Though 
there would be several aspects in the study of the forma¬ 
tion of these rich structures, we believe the self gravity 
is one of the most intrinsic aspects for the structure for¬ 
mation. 

Therefore, in this paper, we would like to study funda¬ 
mental properties of the self-gravitating system (SGS) 
from a view point of statistical mechanics. We specify 
the SGS as a mass distribution which consists of discrete 
mass points mutually interacting through the force of 
Newtonian gravity. The contact collisions of the ingredi¬ 
ents will be neglected. 

The SGS is well known to be unstable and eventually 
collapses. Then it is widely believed that there is no 
thermal equilibrium in SGS. However relevant structures 
would be observed in the intermediate stage of its evo¬ 
lution well before the collapse. Therefore we may have 
a room to construct statistical mechanics of SGS. Actu¬ 
ally there have been many studies dealing with statistical 
mechanics of SGS jlj]. For example, an approach by de 
Vega et al. || considering canonical/grand canonical en¬ 
semble of SGS in order to explain the scaling relation in 
the interstellar medium seems interesting. The instabil¬ 
ity of SGS requires the inevitable introduction of cutoff 
in their work. We would like to circumvent this situa¬ 
tion allowing complex free energy, whose imaginary part 
provides dynamical information of the system. 

We shall construct statistical mechanics of SGS after 
the model of classical electron-gas system on uniform 
ion background. It has the same square-inverse law of 
force but with opposite signature to SGS. Contrary to 
the ordinary belief, the correspondence, i.e. the replace¬ 
ment e 2 —> — Gm 2 , where e is the charge of electron 
and G is the Newton’s constant, is quite suggestive if 


we identify the temperature T as the velocity disper¬ 
sion of SGS. For example, the Debye wavelength A d cor¬ 
responds to the Jeans wavelength A j, and the plasma 
frequency is related to the inverse of the free fall time, 
LOp := —47t Gmn = — 1 /t 2 j. We believe that these cor¬ 
respondences are not merely appearance. Actually, we 
can calculate the free energy of SGS in the same way as 
the classical electron gas system by applying the cluster 
expansion method 11- Moreover, n-point correlation 
functions are similarly obtained. 

b. Cluster expansion method We consider a non- 
relativistic gas system of N particles with the same mass 
m in a box of size L. The gas is in thermal equilibrium 
at T, where the temperature T is identified with the ve¬ 
locity dispersion (v 2 ) of the system, i.e. ksT = m(v 2 )/ 3, 
where ks is the Boltzmann’s constant. 

We work in the canonical ensemble with the Hamilto¬ 
nian of the system: 
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where (f>ij is a potential through which each particles in¬ 
teract. The partition function for fixed number of par¬ 
ticles N, volume V = L 3 , and temperature T, is given 
by 
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where P = (ksT) h is the Planck’s constant, and 
<fiij := -P<j>ij. We have introduced a source function 
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u(xi) for convenience of calculation and integrated over 
the momenta. By using functional derivative with re¬ 
spect to the source u(xi ), fc-point correlation functions 
are given by 


G{x i, • • • ,Xk) = lim 
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For the calculation of thermodynamic quantities, we 
introduce the logarithm of the average of the configura¬ 
tional sum setting u —> 0: 


Wo := W(0) = In 




where (•••}= V N f ■ ■ ■ d N x. 

The cumulant expansion of Wq is possible: 
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where (■ • -) c is a cumulant and 1,2, • • ■ ,M = N(N— 1)/2 
are all the possible pairs of 7V-particles. 

It is convenient to use graphical representations in or¬ 
der to calculate Wo. An each term in Eq.(|>j) corresponds 
to a graph which consists of vertices and lines. Each 
line terminates at two distinct vertices. A vertex Xi sim¬ 
ply represents an integration position and a line connect¬ 
ing Xi and Xj represents an interaction i pij . All the el¬ 
ements in a graph are multiplied with each other and 
integrated over all vertices. This is the cluster expan¬ 
sion method. All the graph which contributes to In Z is 
connected and does not include articulation points which 
divide the graph into plural pieces, because of the char¬ 
acter of the cumulant and the translational invariance of 
the system. 

By expanding with respect to the number density n = 
N/V , Wo is reduced to the following form !!• 
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where Vij is the number of interaction lines between the 
vertices i and j. 

c. Application to the self-gravitating system Let us 
consider a non-relativistic self-gravitating system where 
<j>ij = — Gmf/rij with = \xi — Xj |. 


c.l. Free energy We shall include the higher orders 
in n for calculating Wo- Here we choose a series of graphs 
which is the lowest order of <p in each set of graphs con¬ 
taining k internal vertices. This is equivalent to a sum 
of all graph which has the topology of a ring (ring ap¬ 
proximation shown Fig.[I] (a)). This approximation will 
not be valid for short distances where p grows without 
bound. 
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FIG. 1. A series of cluster expansion: (a) ring graphs, 
(b) chain graphs, (c) multi-lines of the chain graphs. A 
symmetric factor is written in the left side of each graph. 


A ring graph which contains k vertecies in Eq.(jTj) cor¬ 
responds to 
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where p(q) = Anf3Gm 2 /q 2 is the Fourier transform of 
p(x). ( k — l)!/2 terms of this kind for each k yield 


Wn = 


V_ 

7 r 2 

1 


\q\>2n/L 

f-OO 

'2?r/L 


dq ^ [np{q)] h 
( 27r ) 3 k 


= 4^1 dq ">‘ 


K 2 

( k 2 A1 

- 7T — In 

1-— 1 

q 2 

A ( l)\ 


12tt 2 L 

+87t 3 In 


2itk 2 L 2 + k 3 L 3 In 
k 2 L 2 


kL + 27t 


kL — 2tt 


1 - 


47T 2 


+iiT (k 3 L 3 + 87t 3 ) 9 ( kL — 2n) 


(9) 


where k 2 := ATvf3Gm 2 n = 12 t: 2 /\ 2 with A j = 
y/n(v 2 )/(Gmn) being the Jeans wavelength and 9 is the 
Heaviside function. For the integration in Eq.(l|), we an¬ 
alytically continue e 2 to —Gm 2 from the above complex 
e 2 plane by the reason under Eq. ©■ 

We notice that the imaginary part appears for L > 
27 t/k = A j/VS. This originates from the negative ar¬ 
gument in the logarithm for small momentum q. The 
appearance of it apparently indicates the Jeans instabil¬ 
ity. The thermodynamic limit N —> 00 , V —> 00 with 
fixed number density n, yields the free energy: 
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The above choice of analytical continuation guarantees 
the positivity of the imaginary part of F: The system is 
truly dissipative rather than anti-dissipative. 

In general, the imaginary part of a free energy is di¬ 
rectly related with the decay strength of the system T [[| : 
T = (w/3/7r)Im F, where — us 2 is the negative eigenvalue at 
the saddle point which divides a metastable region from 
a stable region. If we identify us as the inverse of the free 
fall time, i.e. us = tJJ = \JAnGmu, the decay strength T 
becomes 
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which is essentially the inverse of the binary relaxation 
time 
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except replacing In N with N. 

Contrary to the extensive variables in the conventional 
thermodynamics, the imaginary part of F apparently is 
superextensive since Imf oc n 3 / 2 but not oc n. More¬ 
over the imaginary part is related to the fluctuation of 
the system through the fluctuation-dissipation theorem. 
In this context, it seems interesting to notice that the 
number 3 in the above n 3 / 2 comes from the spatial di¬ 
mensionality and 2 from the inverse-square law of the 
gravitational force. This reminds us of the Holtsmark 
distribution of the gravitational force acting in the uni¬ 
form self gravitating system |Ml or the stable distribution 
of index 3/2 §. 


c.2. Two point correlation function Let us now 

turn our attention to the correlation functions. One- 
point correlation function is simply a number: G i(5f) = 
N. The normalized two-point correlation function which 
is usually used in astrophysics is 

£(r) = G(x i,f 2 )/ [' G(xi)G(x 2 )], (13) 

where r = \x\ — x 2 \. A similar cluster expansion for 
G(xi,x 2 ) would be 
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We shall include higher orders in n as previous calcula¬ 
tion. Therefore, among each set of graphs which contain 


k internal vertices (excluding both X\ and x 2 ), we choose 
the lowest order skeleton graph in ipij . This is equivalent 
to a sum of all the graph which has the topology of a 
chain shown in Fig.[j]-(b). 
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Since we have analytically continued from the above com¬ 
plex e 2 plane in Eq.(j|), we should choose the pole at 
—nr for the integration. Summing over multi-lines of 
chain graphs shown in Fig{j]-(c), we obtain the following 
function: 

<?i 2 = F \2 + (W12) 2 + ^1 (^12) 3 + “jy (W12) 4 d - 

= exp(^ 12 ) - 1. (16) 


Moreover, we should consider mass renormalization for 
the external vertex. In principle, this is given by sum 
of all the graph which consists of a cluster around the 
external vertex x\. Since infinite summation of this class 
of graphs at present is technically impossible, we phe¬ 
nomenologically introduce an effective mass. Since the 
gravitational attractive force balances with the stirring 
force arising from the velocity dispersion at the length 
scale K -1 , the mass inside of this scale is thought to be¬ 
have collectively. Thus we estimate the effective mass 
m* = (47r/3)(« -1 ) 3 nTO, which should be replaced with 
the mass m at Hi or at x 2 but not both. Strictly speak¬ 
ing, m* is a parameter in the present phenomenologi¬ 
cal argument. By replacing m at X\ in Eq. (pb|) with 
m* at X \, the two-point correlation function becomes 
G{xi,x 2 ) = N 2 (pi 2 \ m • Finally we obtain the 

at x y at x y 

normalized two-point correlation function: 


£(r) = exp 
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The scaling property is manifest: Small scale SGS with 
large n has the same correlation function as that of large 
SGS with small k. 

This function £(r) has interesting features. As a func¬ 
tion of s = nr, it has a unique inflection point at 
s = 0.440 when plotted in the Log-Log graph. The slope 
at s = 0.44 is —1.667 and the magnitude there is 0.986. 
Therefore, the inflection point r c ~ 0.44k -1 ~ 0.04Aj is 
regarded as the correlation length of SGS. 
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FIG. 2. A Log-Log graph of two-pint correlation func¬ 
tion: £(s) = exp [cos(s)/(3s)] — 1 (solid line), (s/0.45) -1 ' 77 
(broken line), observational data of galaxies H] (open circle), 
and observational data of clusters of galaxies [IL(| (filled cir¬ 
cle). These correlation lengths of the observational data are 
rescaled. 


above values through two-point correlation function we 
obtained. 

d. summary and outlook Studying the canonical 
ensemble of the self-gravitating system (SGS), we ob¬ 
tained the complex free energy by summing an infinite 
series of graph in the cluster expansion method for SGS, 
Eq. ©• The imaginary part of the free energy yields the 
decay strength of SGS, Eq. ©■ Similar summation of 
an infinite series of graph yields the universal two-point 
correlation function Eq. (|l7|) which scales essentially with 
the Jeans wavelength. The correlation length is linearly 
proportional to the mean separation of ingredients. 

We hope further development will be reported soon, 
including (a) the higher-point correlation functions, (b) 
much profound calculation on the complex free energy, 

(c) systematic argument on the mass renormalization m*, 

(d) the observational tests of our arguments, and (e) the 
effects of the cosmic expansion and of the Dark matter. 


The above two-point correlation function is shown in 
Fig.|| with typical observational results {~||TT| on top of it. 
The observation does not exclude our correlation function 
except small scale region where the interaction ip exceeds 
1 and the approximation we used is obviously invalid. 

From the observational data the correlation 

lengths of galaxies and of clusters of galaxies directly read 
off as r ya,axy ~ 6.2/i -1 Mpc and r f ueter ~ 15/i -1 Mpc re¬ 
spectively. Since r c ~ 0.44k -1 for the two-point correla¬ 
tion function we obtained, we rescaled the observational 
data: r —> s = nr with k = 1/14 /iMpc -1 for galax¬ 
ies and k = 1/34 ft.Mpc -1 for clusters of galaxies, where 
h = Hq/ (l00knrs -1 Mpc -1 ) with Hq being the Hubble 
constant at the present time. These values of k corre¬ 
spond to ~ 152/i -1 Mpc and A“ i “ ster ~ 370/i -1 Mpc 

from the relation given immediately after Eq.([|). On 
the other hand, typical Jeans lengths calculating from 
the standard observations for galaxies and for clusters of 
galaxies are respectively 


\ a ; laxy = 123Mpc 


( W )( 

Y 500kms 1 J V 


3.7 x 10 -2 /i 3 Mpc 


2 x 10 44 g 


m 

-3 


x 

1/2 


fluster = 450 M P C 


(, 2 ) 


lOOOkms” 


2 x 10 47 g 


1.1 x 10 -5 /i 3 Mpc 


-3 


- 1/2 


(18) 


(19) 


These values of the Jeans lengths are not so far from the 
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